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ABSTRACT 

The cold dark matter (CDM) scenario generically predicts the existence of triaxial dark matter 
haloes which contain notable amounts of substructure. However, analytical halo models with 
smooth, spherically symmetric density profiles are routinely adopted in the modelling of light 
propagation effects through such objects. In this paper, we address the biases introduced by 
this procedure by comparing the surface mass densities of actual N-body haloes against the 
widely used analytical model suggested by Navarro, Frenk and White (1996) (NFW). We 
conduct our analysis in the redshift range of 0.0 - 1.5. 

In cluster sized haloes, we find that triaxiality can cause scatter in the surface mass density 
of the haloes up to (7+ = +60% and ct- = —70%, where the 1-cr limits are relative to the 
analytical NFW model given value. Subhaloes can increase this scatter to (t+ = +70% and 
(T_ = —80%. In galaxy sized haloes, the triaxial scatter can be as high as (7+ = +80% and 
(T_ = —70%, and with subhaloes the values can change to (7+ = +40% and (t_ — —80%. 

We present an analytical model for the surface mass density scatter as a function of dis- 
tance to the halo centre, halo redshift and halo mass. The analytical description enables one 
to investigate the reliability of results obtained with simplified halo models. Additionally, it 
provides the means to add simulated surface density scatter to analytical density profiles. As 
an example, we discuss the impact of our results on the calculation of microlensing optical 
depths for MACHOs in CDM haloes. 
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1 INTRODUCTION 

The cold dark matter model, in which the non-baryonic part of 
the dark matter is assumed to consist of particles that were non- 
relativistic already at the time of decoupling, and that interact pre- 
dominantly through gravity, has been very successful in explain- 
ing the formati on of large-scale structures in the Universe (see e.g. 
|Primacldl2003l , for a review). In this scenario, both galaxies and 
galaxy clusters are hosted by CDM haloes, which formed hierar- 
chically through mergers of smaller subunits. 

Even though N-body s imulations generic ally predict CDM 
haloes to be triaxial (e.g. J ing & Sut with substantial 

amounts of substru ctures left over from the merging process (e.g. 
lMooreetal.lll999l) . simplified halo models are often adopted in 
the modelling of light propagation through such objects. The most 
common approach is to treat dark matter haloes as spherical ob- 
jects with smooth density profiles, usually either of the NFW 



I Navarro. Frenk & Whitelll996h form, some generalization thereof 
i Zhadll996h . or that of a cored or singular isothermal sphere. 

The light emitted from high-redshift objects such as quasars, 
supernovae, gamma-ray bursts, galaxies and galaxy clusters will 
typically have to pass through many dark matter haloes before 
reaching an observer on Earth. Several investigations have al- 
ready indicated that smooth and/or spherical halo models may 
lead to incorrect results when treating the gravitational lens- 
ing effects associated with such foreground mass condensations 
(e.g iB artelmann & Weiss' 'l994'; 'Dalai. Holder & Hennawil |2004| : 
lOguri & Keeton 2004 ; .Hennawi et al. 2007) 

More realistic features like triaxiality and substructures can 
be included in gravitational lens calculations either by employ- 
ing N-body simulati o ns directly (e.g. B artelmann & Weiss 19941: 
ISeliak&Hol3ll999l : iHolopainen etal] 1200^ or by using ana- 
lytical expressions f or the halo sha pes (e.g. iKochanelj ll993 : 
iGolse & Kneibll2002l : [Evans & Hunt ei1l200d: Bia e 2003) andsub- 
halo properties (e.g. lOguril 120051 : IZackrisson & Riehml l2007h . 
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While N-body simulations often represent the safest choice, the 
approach is computationally demanding and does not always al- 
low one to identify the features of the mass distribution responsible 
for a specific lensing effect. Methods which bring simple, analyt- 
ical halo models into contact with the full phenomenology of the 
N-body simulations are therefore highly desirable. 

In this paper, we focus on the projected mass density of 
CDM haloes as a function of distance from the halo centre. There 
are several situations in gravitational lensing when realistic es- 
timates of the surface mass density (i.e. convergence) along a 
given line of sight through a dark halo may be important. Ex- 
amples include the calculation of image separations in strong 
lensi ng by subh aloes located in the external potential of its host 
halo l lOgurill20b5i' ). attempts to correct the luminosities of super- 
novae type la fo r the magnification by foreground haloes (e.g. 
lGunn arsson''2004) and estimates of the distribution of microlensing 
optica l depth s for high-redshift MACHOs (e.g. IWvithe & Turned 
I2OO2I ; IZackrisson & Riehdl2007h . Other applications include the 
assessments of light propagation effects in models with non- 
zero coupling betweeen da rk matter particles and photons (e.g. 
IProfumo & Sigurdsonll2007h . 

Here, we use high-resolution, dissipationless N-body simula- 
tions of CDM haloes to investigate the errors in surface mass den- 
sity introduced by treating these objects as spherical with smooth 
density profiles of the NFW type. Simple relations for the surface 
mass density error as a function of halo redshift and distance to the 
halo centre are presented, making it easy to investigate the reliabli- 
tiy of results obtained with simplified halo models . 

On a related note, iKnebe & Wiessne3 ( l2006h recently inves- 
tigated the error introduced by spherically averaging an elliptical 
mass distribution. They found that for axis ratios typical for cosmo- 
logical dark matter haloes, the variance in the local density can be 
as large as 50% in the outer parts. The current paper examines the 
problem of halo triaxiality from a slightly different point of view. 

The N-body simulations used are described in Section 2. In 
Section 3, we describe the methods for extracting the halo sample. 
In Section 4, we compare the CDM surface mass densities obtained 
along random sightlines through the N-body haloes to the corre- 
sponding results obtained from smooth and spherical NFW models 
fitted to the same haloes. Section 5 presents a set of simple relations 
for the surface mass density errors introduced by this procedure as 
a function of distance to the halo centre and halo redshift. Section 6 
discusses how these relations may be used in the context of optical 
depth estimates for MACHO microlensing. A number of caveats 
are discussed in Section 7. Section 8 summarizes our findings. 



2 N-BODY SIMULATIONS 



2.1 The 64 h ^Mpc simulations 

For studying cluster sized haloes, we utilize a suite of four high- 
resolution N-body simulations. The simulations were carried out 
using the publicly available adap tive mesh refinement code MLAPM 
jKnebe. Green & BinnevI 1200 ll ). focusing on the formation and 
evolution of dark matter galaxy clusters containing of order one 
million particles, with mass resolution 1.6 x 10** M0 and 
spatial force resolution ~2/i~^ kpc. They are so-called "zoom" or 
multimass simulations in which we first created a set of four inde- 
pendent initial conditions at redshift z = 45 in a standard ACDM 
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Figure 1. The masses and redshifts of the haloes used in our analysis. The 
total number of haloes is 336. See text for details. 

cosmology (fio = 0.3, = 0.7, fif, = 0.04, h = 0.7, erg = 0.9). 
512'^ particles were placed in a box of side length 64h~^ Mpc giv- 
ing a mass resolution of nip — 1.6 x W^h^^ M©. For each of 
these initial conditions we iteratively collapsed eight adjacent par- 
ticles to a single particle, reducing our simulation to 128"^ particles. 
These lower mass resolution initial conditions were then evolved 
until z — 0. At z = 0, eight clusters from different regions of our 
simulations were selected: 4 halos from box #1, one from box #2, 
one from box #3 and 2 from box #4. The masses of these haloes 
are in the range I-3x 10^*^"^ M© and triaxiality param eters vary 
from 0.1 to 0.9. Then, as described by iTormeg il997h . for each 
cluster the particles within five times the virial radius were tracked 
back to their Lagrangian positions at the initial redshift (z — 45). 
Those particles were then regenerated to their original mass resolu- 
tion and positions, with the next layer of surrounding large particles 
regenerated only to one level (i.e. 8 times the original mass resolu- 
tion), and the remaining particles were left 64 times more massive 
than the particles resident with the host cluster. This conservative 
criterion was selected in order to minimise contamination of the 
final high-resolution haloes with massive particles. 

A more elaborate description of this data set and a de- 
tailed investigation of the sense of rotation of the satellites and 
the properties of the tidally induced debris field of disrupt- 
ing satellites can be found elsewhere jWarnick & Knebd l2006l : 
IWarnick. Knebe & Powe3l2007 dlbl) . 

2.2 The 10 /i ^ Mpc and 40 Mpc simulations 

For studying galaxy sized halos (A/ ~ 10^"^ Mq) and for acquiring 
better statistics on the larger haloes, we ran two additional simula- 
tions with smaller box sizes. The same simulation code MLAPM as 
in the 64 h~^Mpc simulations was used, but the cosmological con- 
stants were slightly different: Qo = 0.27, Qb = 0.044 and fix = 
0.73. These two simulations are standard cosmological simulations 
where all the particles have the same mass: 4.47 X lO^/i"^ Mo 
and 2.86 x lO^h'^ M© for 10 h'^Mpc and 40 h'^Mpc simula- 
tions, respectively. Spatial force resolutions were 0.46 h,~^kpc and 
1.8 ft^^kpc. Both simulations were followed until 2 = and halos 
were identified at different redshifts. The 10 ft^^Mpc simulation 
was followed from z — 71.52 and the 40 h^^Mpc simulation from 
z = 47.96. 

One may argue that the large scale modes, ignored in these 
relatively small volume simulations, can cause spurious errors in 
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the results. Especially in the smallest simulation volume, the long 
wavelength perturbations will not be present. However, the distri- 
butions of mass concentrations remain the same when compared 
to larger volume simulations. Also, in the 10 /i^^Mpc box, we re- 
strict our analysis to intermediate mass halos and their subhalos. 
We believe that the simulations are reliable for our purposes and in 
agreement with simulations covering larger volumes. 



3 CDM HALOES 

3.1 Finding and truncating haloes 

Finding and truncating dark matter haloes within cosmological 
simulations is an interesting and challenging task. Many authors 
have addressed this problem successfully by developing sophisti- 
cated algorithms which can locate hal oes by a variety of techniques 
(e.g. lDavisetalJll985l: iFrenk et al] 1 19881: iBertschinger & Gelb 
1991; Suto, Cen & Ostriker 19921 ; IWeinberg. Hernquist & Katz 
1997.; Klypin & Holtzman 199^7^ 

Our analysis examines the cores as well as the outskirts of the 
haloes without the luxury of being allowed to overlook the exact 
properties of the low density regions. We need to study all parts 
or the halo density profile in three dimensions, and we need to get 
rid of the background particles at all radii as well as possible, espe- 
cially near the virial radius. To achieve this, we use a highly capable 
halo finder which can determine the parent potentials of individual 
particles and "clean" host haloes from their subhaloes. 

We find and truncate our haloes and su bhaloes using the 
MLAPM Halo Finder (MHF) dOill et alj l2004t) . MHF uses the 
adaptive grids of MLAPM to locate haloes within the simulation. 
MLAPM's adaptive refinement meshes follow the density distribu- 
tion by construction. The grid structure naturally "surrounds" the 
haloes, as the haloes are simply manifestations of over-densities in 
the mass distribution of the simulation box. The grids of MLAPM are 
adaptive, and it constructs a series of embedded grids, the higher 
refinement grids being subsets of grids on lower refinement levels. 
MHF takes this hierarchy of nested isolated grids and constructs 
a "grid tree". Within that tree, each branch represents a halo, thus 
identifying haloes, subhaloes, subsubhaloes and so on. 

While a branch of the tree identifies the majority of particles 
associated with a halo, the surrounding region is checked for addi- 
tional particles if the halo is embedded within another halo (i.e., it is 
a subhalo). To gather additional particles, a larger collection radius 
is defined, and all the particles within this radius are assigned to the 
halo. In this paper, the collection radius has been defined as half the 
distance between the current halo and the next most massive halo. 

The gravitationally unbound particles are then removed from 
the haloes in an iterative fashion. If a particle is not bound, it is 
assigned to the subhalo's host or the background as appropriate. 
This, however, does not guarantee that each particle is uniquely 
assigned to a halo. It is possible for a particle to be shared by two 
or more halo potentials by these criteria. 

Using this set of particles, the canonical properties of the 
haloes are calculated. For example the virial radius is found by 
stepping out in (logarithmically spaced) radial bins until the density 
reaches phaio(7'vir) = Avir(z)p6(z), where pt is the cosmic matter 
density. Particles outside this radius are removed. If this density is 
not reached, then we consider the furthest bound particle from the 
centre of the halo as the radius. 



Table 1. Summary of the number of haloes and the particle counts in our 
sample. First column gives the simulation box size and A'^halocs refers to 
the number of haloes the sample contains per redshift. N^^^ refers to the 
number of pai'ticles within Tvir of the least massive halo and N^]^'^ to the 
number of particles within r^i^ of the most massive halo, in each box. 
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3.2 The halo sample 

Typically, MHF finds thousands of haloes within a simulation box, 
but the number of haloes which have sufficient mass resolution for 
our purposes is unfortunately quite low. When choosing the halo 
sample, we have to balance between particle resolution and the 
number of haloes. After looking at the consistency of the NFW 
fit and the halo mass profiles, we decided to include the ten most 
massive haloes from the 10 /i~^Mpc and 40 h'^Mpc simulations. 
This choice translates to a minimum resolution of 36,000 particles 
in the smallest halo in our sample. Note that our particle counts are 
given without subhaloes, which normally contribute about 10 — 20 
% of the host's mass. 

Nearly all haloes in the sample contain more than 100,000 par- 
ticles at the maximum redshift, Zmax = 1-5. The only exception is 
the 40 h'^Mpc box, in which the number of particles contained by 
the tenth most massive halo is only 36,000. However, this is com- 
pensated by the large particle counts in the multimass simulations, 
which contain haloes in the same mass range. The particle counts 
of the halo sample is shown in Table[T] Note that the largest number 
of particles within a halo is always found at 2; = and the smallest 
number of particles at z = 1.5. 

Each of the eight 64 /i~^Mpc multimass simulations contains 
only one high resolution cluster halo region. When the virial mass 
accretion histories were examined, we found that the youngest of 
the eight clusters is composed of three smaller haloes with compa- 
rable masses. In this case, it was impossible for MHF to separate 
the three interacting haloes from each other, and therefore, we were 
left with seven well defined haloes per redshift from the multimass 
simulations. Excluding the "triple-cluster" does not introduce any 
selection bias since this type of specific merger occurs rarely. 

The number of redshifts used per simulation box can be seen 
in Fig.[T] The redshift counts are 13, 15 and 8 for the 10 /i~^Mpc, 
40 h^^Mpc and 64 h^^Mpc boxes, respectively. In total, the num- 
ber of haloes in the sample is 130 -1- 150 -1- 56 = 336. 

Our halo sample is divided into two mass classes with a mass 
gap at M ~ 10^^ MqH'^ at 2; = (see Fig.[T]l. We refer to the 
more massive class as the cluster sized haloes and to the lower mass 
class as the galaxy sized haloes. We would need another cosmolog- 
ical simulation with a comovins box size of ~ 25 h-^Mpc to fill 
the mass gap. The 40 h~^Mpc simulation does contain haloes with 
masses at the range of the gap, but unfortunately the resolutions of 
these haloes are insufficient for our purposes. 
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3.3 Analytical haloes 

Once we have found, extracted and chosen our N-body haloes, we 
want to construct their analytical counterparts. We have chosen to 
use the NFW density profile so that our fit procedure can be com- 
pared to earlier work and reproduced easily. The fit equation is 

ao 



p{r) = 



-(1 



(1) 



where ao and ai are our fit parameters. We measure the number 
density of the particles in a halo at a certain radius by dividing 
the particles in logarithmically spaced radial bins (shells) and cal- 
culating the particles within these bins. Then, by the fit procedure 
(Levenberg-Marquardt algorithm), we find the ao and ai values 
which minimize 

2 



X 



(2) 



where Ps{ri) is the measured number density of the ith shell and 
ri is the midpoint radius of that shell. The estimated "measurement 
error" in the number of particles within a shell is assumed to be 
Poissonian: ~ VNi- This makes the estimated "measurement 
error" in the number density to be (jp(ri) — Ni/ p^{ri). The choice 
of weighting the fit with o-p(ri)^ tends to provide the inner core less 
weight than to the outer regions. This is a deliberate choice because 
our analysis is affected by the full structure of the halo and only in 
few rare cases dominated by the center. 



4 MEASURING THE SURFACE MASS DENSITY 
4.1 Sightlines 

After we have assigned the particles to a certain halo, we smooth 
the mas s distribution by treating the p articles as triangular-shaped 
clouds jHocknev & Eastwood 1 19811) . This approach allows us 
to measure the surface mass densities of the haloes using one- 
dimensional sightlines. If the particles were treated as point masses, 
the surface mass density would have to be measured using tubes of 
finite radius, which would lead to an average measurement within 
a tube. To avoid the averaging process, each particle is assigned a 
surface mass density 



Sp(r) 



Eo(l- 




r < a 
r > a 



(3) 



where r is the impact parameter between a sightline and a particle 
and A/p is the mass of the particle, a is the size of the adaptive grid 
cell where the particle is found, and a is used in our analysis as the 
effective radius of the particle. The cell size depends from the local 
particle density and is determined by the cosmological simulation 
grid constructing process (MHF and MLAPM can both construct the 
same adaptive grid tree). The cell sizes in the simulations can range 
from 250 /i~^kpc to 0.150 /i~^kpc, and the smallest cells are au- 
tomatically located in the dense cores of haloes and subhaloes. 

To sample the surface mass density of an MHF extracted halo, 
we pierce it with 10,000 randomly oriented sightlines with impact 
parameters chosen from a uniform random distribution in the range 



G (0, 



We then calculate the surface mass density for 



each sightline EsL(7'imp) by summing the surface mass densities 
of the particles for which r < a. 

The proper truncation of the smooth mass distribution at the 
virial radius is slightly complicated. Some of the particles can ex- 
tend a portion of their mass beyond the virial radius. For these 



particles, only the mass inside r^ir is taken into account. Further- 
more, the extended particles introduce an additional requirement 
for our halo finder - it needs to be able to find particles which 
belong to a given halo potential (based on the particle velocity) 
out to ~ 1.5rvir. This is because most of the particles between 
1.0 < r-vir < 1-5 extend their effective radii (defined in Eq. [3]l 
into the region inside the virial radius and contribute mass to our 
sightlines. 



4.2 Truncated density profile 

We want to analyse how the surface mass density of an N-body halo 
behaves compared to an analytical model. To calculate the analyti- 
cal value, we use a truncated version of the projected NFW density 
profile, which gives us the predicted surface mass density as a func- 
tion of the impact parameter (rimp) and the virial radius of a halo 
(r-vir). Truncation is needed because our haloes have a limited size, 
unlike the NFW haloes, which extend to infinity. The fit parameters 
ao and ai are also used in the truncated model - after they have 
been determined by fitting Eq.[T]to the respective halo. 

We derive the following form for the truncated NFW surface 
mass density 



SNFw(a;,c)= / p{z)dz = 2aoaiF{x,c), 



(4) 



where 



F{x,c) 



Here x 



-1 c+x^ 
cosh ( I . ) 
^ (l + c)3: ' 



1 + C 



1 + C 




X < 1 



X > 1. 



(5) 



p/ai and c 



unc/ai. Note that c is analogous 



to the NFW concentration parameter Cvir when the truncation ra- 
dius is equal to the virial radius rvir. When |a; — 1| < 0.1, we 
have to interpolate (linearly) between both forms of the function, 
i^x>i(l.la;, c) and Fi<i(0.9a;, c) for numerical reasons. 

For the untruncat ed version of Eq.|4l in which c ^ oa, see e.g. 
iGolse & KneibI ( l2002h . The untruncated model extends asymptot- 
ically to infinity and does not have the turnover seen for example 
in Figure [2] between SOOkpc < r < 1300kpc. Th e surface mass 
density has also been derived bv lBartelman"ril ( ll996l) . 



5 A SIMPLE CORRECTION SCHEME FOR SMOOTH 
AND SPHERICAL DARK MATTER HALOES 

5.1 Fitting the scatter 

Figure 2 shows the surface mass density scatter of a typical cluster, 
with its subhaloes removed, and the scatter is seen to follow the 
NFW profile quite well. For some haloes, there are some strong 
deviations within the inner 5 % of the virial radius, and for the 
majority, a weak, continuous offset trend near the virial radius is 
seen. The problems o f the NFW profile in the core is a known issue 
( lNavarroetZll2004l) . The scatter in surface mass density around 
the mean is roughly log-gaussian, as will be seen. 

The scatter seen in Figure [2] is a consequence of the fact that 
the haloes are not spherically symmetric ~ if they were, all the 
measured points would land on the NFW profile curve. Subhaloes, 
when not excluded as in Figure (2] cause a small number of sight- 
lines to produce even higher surface mass density values than the 
ones seen in Figure|2l 
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Figure 2. This figure shows how the surface mass density scatters due to 
triaxiaUty in most of the haloes in our halo sample. Here the values are mea- 
sured from a cluster sized halo with 10,000 randomly oriented sightlines, 
without subhaloes. The dashed line shows the truncated NFW fit profile, and 
the two vertical lines mark the borders of a radial bin. The SgL-distribution 
of this bin is fitted with a log-normal function in Fig. [5] 




Figure 3. A typical log-normal fit of the scatter. The sightlines have been 
taken from the halo shown in Fig. [2] from an impact radius range of 0.7 < 
^imp/^vir < 0-8- The Surface mass densities are divided in 30 bins, of 
which only 1 8 bins containing the majority of the sightlines are visible here. 
The log-normal distribution is fitted to the frequency of the sightlines in 
these 30 bins. All fits in our analysis are done accordingly. 



For comparing the surface mass densities on sightlines from 
different haloes (with different virial radii and masses), the impact 



parameter is transformed to a unitless variable x\n 



The amount by which each EsL(a;imp) deviates from the 
NFW profile measures how much the shape of the halo deviates 
from a spherically symmetric analytical model at a specific sight- 
line. This is why the surface mass density is transformed relative to 
the NFW profile of each halo as follows: 

SsL(a::imp) 



'S'{3^imp) — 



(6) 



SNFw(3;impCvir, Cvir) ' 

where Cvir = rvir/ai and Enfw(k, c) is the function in Eq.|4] As 
a result of all this, the impact parameter has been scaled to the units 
of the virial radius and the measured surface mass density to the 
units of the NFW value. 

We divide the measurements in 10 radial bins (one example 



bin is shown by the vertical lines in Fig.|2}. The width of each bin 
is Aa:imp = 0.1 and each bin contains 1000 measurements. These 
radial bins are further divided in 30 bins in surface density, which 
typically contain a maximum of ~ 200 measurements in a single 
bin (see Fig.O. 

Once the measurements points are binned, we fit log-normal 
distributions to the S^Xiynp) distribution, for each radial bin, indi- 
vidually for each halo. The fit function is 

-(ln{S)-M)V{2<T^) 



(V) 



Repeating the fit procedure for all haloes in the sample allows 
us to derive the best fit parameters as a function of impact parame- 
ter and the halo redshift: o-(a;iinp, z) and ^{ximp, z). Essentially, cr 
measures the width of the scatter and measures the mean devia- 
tion from the NFW fit profile in ln(S')-space. Because each halo is 
divided in 10 impact parameter bins, the total number of log-normal 
fits (and both fit parameters) is 336 x 10 = 3360. 

5.2 Finding dependencies to z and ximp in the log-normal 
parameters 

Finally, we searched for trends in the ^(simp, 2, A/vir) and 
/i(a;imp, z, Mvir) data for constructing an analytical description of 
our measurements. As the analysis progressed, we quickly became 
aware of the fact that the log-normal parameters do not correlate 
with the virial mass of the halo as strongly as with the other two 
variables. Thus, we reduced our analytical description to functions 
o"(a;imp, z) and /^(ximp, z) for the galaxy and cluster sized haloes 
separately. We also give the analytical description with and with- 
out subhaloes, for both mass classes, which then makes the number 
of our final analytical descriptions four. All four descriptions are 
expressed by changing parameters within the following function 
forms: 

(T(a;imp, z) = Po(2;imp) + zPl{Xiuip), (8) 
^i{ximp,z) = Qo{xin-ip) + ^Ql (aJimp ) , (9) 

where both Pi(a^imp) and Qi(ximp) are second order polynomials: 



[Xi^ 



— <?i2a;imp + IJilKimp + (?iO • 



(10) 
(11) 



Examples of the analytical descriptions of ^(ximp, z) and 
n{xijn-p,z) with the actual data can be seen in Figures |4] and [5] 
The constants for all polynomials are given in Tables |2] and |3] Be- 
fore the analytical versions of CT(a::iinp, 2) and /^(ximp, z) are fitted 
to the log-normal parameters, the data are averaged over bins with 
a width of Az ~ 0.2 to reduce noise. Also, we had to drop the 
inner 5 % of the data (a;imp < 0.05) because within this region, the 
NFW fit fails to follow the data correctly in a significant number 
of cases. This is because of our choice of fit weighting (see Section 
[33J. 

Equations[8]and|9]quantify the distribution of the surface mass 
density of a generalized galaxy or cluster sized CDM halo at any 
impact radius and redshift within the ranges we have used. Because 
of the log-normal distribution, the geometric mean and standard 
deviations at a given a::imp and z values in ^(ximp, ; 
oM(a;imp,z) 



^J.s{ximp, z) 



, z )-space are 

(12) 

o-s(a:in.p,2) =e"("™P'"\ (13) 
and the upper and lower 1-(t limits in Sixi-^^, 2)-space are 



6 J. Holopainen et al. 



Table 2. The constants for polynomials Pq and Pi , which make up the 
analytical description for cr{x\^p, z) (Eq.[8). The samples are coded as fol- 
lows. G-ns: Galaxy sized haloes, no subhaloes. G-ws: Galaxy sized haloes, 
with subhaloes. C-ns: Cluster sized haloes, no subhaloes. C-ws: Cluster 
sized haloes, with subhaloes. 



Sample 


POO 


POl 


P02 




PlO 


Pll 


P22 


G-ns 


0.217 


-0.0692 


0.305 




0.0514 


-0.0889 


0.289 


G-ws 


0.211 


-0.0659 


0.380 




0.0424 


0.199 


-0.00372 


C-ns 


0.272 


0.0355 


0.299 




0.0386 


-0.121 


0.269 


C-ws 


0.287 


0.0861 


0.266 




0.0260 


-0.00533 


0.247 



Table 3. The constants for polynomials Qq and Qi, which make up the 
analytical description for /^(xiinp, 2) (Eq.|9j. The samples are coded as in 
TableU] 



Sample 


900 


goi 


902 




qio 


qii 


922 


G-ns 


-0.0322 


0.286 


-0.526 




-0.0739 


0.0406 


0.0376 


G-ws 


-0.0612 


0.254 


-0.530 




-0.060 


0.0365 


-0.135 


C-ns 


-0.0333 


0.148 


-0.433 




-0.0746 


-0.00968 


0.0296 


C-ws 


-0.0739 


0.218 


-0.593 




-0.0198 


-0.144 


0.122 



atixi^p, z) = e(M(-.„p,^)+'^(-.„p,^)) (14) 
aj(xi^p,z) = e(''(-'-P'-)-(-™P'-». (15) 

We strongly recommend that the analytical function form 
which includes subhaloes is used with care. This is because sub- 
haloes can have unwanted effects to the log-normal fitting proce- 
dure. The preferred way of using the derived analytical description 
is to use the version which does not include subhaloes and then 
add the subhalo scatter afterwards if really necessary (the statisti- 
cal subhalo contribution to surface mass density is small). A more 
reliable estimate of the subhalo contribution can be acquired by us- 
ing the known subhalo mass and distribution f unctions, which have 
been studied in detail by e.g. lGao et al.l j2004h . 



5.3 Comparing the model to the data 

We confirm the reliability of our fit procedures, the choice of fit 
equations and our analytical description of the data by comparing 
the EsL(a;imp)-values to the predictions of the model. This is done 
for all four models by measuring the frequency of the following 
value: 



0.8 - 




0.2 0.4 0.6 0.8 1 1.2 1.4 



z 

Figure 4. An example of the best fit for the averaged a{x\^p, z) fit pa- 
rameters. Included in this fit are the galaxy sized haloes without subhaloes. 
Each polyline represents data from a single radial bin and is fitted with the 
linear function in Equation [s] rj{xi-^-p, z) increases with growing impact 
parameter at a given redshift. The highest line is the fit for the most distant 
radial bin Simp = 0.95, where 0.9 < I'imp < 1-0. 



0.05 M 




-0.35 ^ ' -^ 

0.2 0.4 0.6 0.8 1 1.2 1.4 

z 

Figure 5. An example of the best fit for the averaged ^l(xi^p, z) fit param- 
eters. The halo sample is the same as in Figure|4] Each polyline represents 
data from a single radial bin and is fitted with the hnear function in Equation 
[9] The absolute value of /^(xiinp, z) increases with growing impact param- 
eter at a given redshift. The lowest line is the fit for the most distant radial 
bin Ximp = 0.95, where 0.9 < aiimp < 1-0. 



r(:Ci^„ z) = ln(g(^imp))-^(^i.np,^)_ ^jg^ 
/^(^imp 7 Z / 

The probability density of T(xiinp,z) should follow a stan- 
dard normal distribution for all redshifts if the model is reliable. 
As shown by the example in Figure |6l this is the case when all the 
10,000 X 130 sightlines for the galaxy sized haloes (without sub- 
haloes) are considered at all redshifts. 

The fact that the resulting distributions are centered on zero 
tells us that /i(a;inip, 2) is reconstructed correctly. The standard de- 
viations of the distributions are close to unity, which again tells us 
that cr(a;imp, z) is a fair estimate of the behaviour of the data at 
all impact radii and redshifts. The test was equally successful for 
galaxy and cluster sized haloes, with and without subhaloes as the 
test shown in Figure |6] 
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Table 4. The minimum and maximum geometric standard deviations and 
mean offsets from NFW in ^(xinjp, z)-space (see Equations 1121 and I13t . 
The samples are coded as in Table[2] For example, for galaxy sized haloes, 
with subhaloes, at z = 1.5 and with impact parameter x^^-p = 1.0, 
the 1-cr upper limit in surface mass density in units of the NFW model 
value is 6^+°" = 0.56 x 2.41 = 1.35, and the lower limit is et^^" = 
0.56/2.41 = 0.23. Note that the S{xi^p, z) distribution is not symmetric 
but log-normal. 



Function 


G-ns 


G-ws 


C-ns 


C-ws 


0-3(0.0,0.0) 


1.24 


1.23 


1.32 


1.33 


0-3(0.0, 1.5) 


1.34 


1.32 


1.39 


1.38 


os(1.0,0.0) 


1.57 


1.69 


1.83 


1.89 


os(1.0,1.5) 


2.29 


2.41 


2.42 


2.82 


Ms(0.0,0.0) 


0.97 


0.95 


0.97 


0.93 


Ms(0.0,1.5) 


0.87 


0.86 


0.87 


0.90 


Ms(l-0,0.0) 


0.76 


0.72 


0.73 


0.64 


Ms(l-0,1.5) 


0.77 


0.56 


0.67 


0.60 




-10 12 3 

T = [ln(i: / Infw) - V'i^' z)] / oC", z) 

Figure 6. Here one of our analytical models (galaxy sized haloes, no sub- 
haloes) is compared against the data. Different curves represent different 
redshifts, and there is no trend between them. All curves are close to the 
standard normal distribution (which peaks at ~ 0.4), which means that our 
model is able to describe the distribution of the measured values quite well. 
See SectionfOlfor details. 



6 IMPACT ON THE MACHO OPTICAL DEPTH 
CALCULATIONS 

MACHOs (Massive Astrophysical Compact Halo Objects) repre- 
sent one class of dark matter candidates, which may be detected 
through gravitational microlensing effects as they pass through 
the line of sight to distant light sources. Although the MACHO 
acronym was originally invented with baryonic objects like faint 
stars and stellar remnants in mind, several non-baryonic dark 
matter candidates l ike axion aggregates (lMembrado|[l998h . mir- 
ror m atter objects jMohapatra & Teplitj|l999[). primordia l black 
holes jGreenl200(]h, quark nugg ets jChandra & Govall200(]h . preon 
stars (Hansson & Sandin 2005) and scalar dark matter miniclus- 
ters ( Zurek, Hogan & Quinn 2006) can also manifest themselves in 
this way. While non-baryonic MACHOs could in principle account 
for a substantial fraction of the dark matter, microlensing searches 
based on observations of light sources in the loc al Universe suggest 
a MACHO dark matter fraction of ~ 20 % (e.g. lAlcock et alj2OO0l : 



ICalchi Novati et alj2005]) . Low-redshift microlensing observations 
are only able to detect the densest, most compact MACHOs at a 
given mass. This does not necessarily represent a robust upper limit 
on the relative importance of MACHOs, and high-redshift observa- 
tions will be required to settle the issue (see lZackrisson & Riehml 
I2OO7I for a more detailed review). 

While difficult to directly relate to observational quantities, the 
microlensing optical depth r is often used for estimates of the via- 
bility and relevance of different microlensing scenarios. Formally, 
the microlensing optical depth represents the average number of 
lenses along a random line of sight. Under the assumption that the 
lenses do not overlap along the sightline, it also represents the frac- 
tion of sky that is covered by regions in which a point source will 
be microlensed. In the limit of small r the optical depth can there- 
fore directly be used as an estimate of the microlensing probability. 
At higher r, this interpretation does however break down because 
of overlapping microlenses. Here, the optical depth will be used to 
discuss the impact of surface mass density variations due to clumpy 
and non-spherical CDM haloes on MACHO microlensing calcula- 
tions. In the following, we assume the spatial distribution of MA- 
CHOs to follow that of the CDM, as would be expected if they are 
non-baryonic. We caution that these results do not necessarily apply 
to baryonic MACHOs, since such objects may have a spatial dis- 
tribution that is substantially different from that of the overall dark 
matte r (as illustrated by population III stars; e.g. lScannapieco et al.l 
l2006h . 

In the extreme case of having all CDM in the form of MA- 
CHOs, the MACHO optical depth along a sightline passing through 
a single halo can be approximated by: 

r=^, (17) 

where Eqdm is the CDM surface mass density of this particular 
sightline and Ec is the critical surface mass density for lensing. 
The latter is given by: 



4^GA.O„i' ^^^^ 

where Doi, Dis and Dos are the angular-size distances from ob- 
server to lens, lens to source and observer to source, respectively. In 
the case when only a fraction of the CDM is the form of MACHOs, 
and the rest is in the form of a smoot h component, the optica l depth 
can instead be estimated using (e.g.' Wvithe & Tumeill2002l where 
an analogous case with stellar microlenses in a smooth matter dis- 
tribution is considered): 

^ ./macho ScDM ^-^gy 

|Sc — (1 — /macho)Scdm| ' 

where /macho represents the MACHO fraction. 

In situations where the sightline is dominated by a single halo 
of (approximately) known mass, and the impact parameter of this 
sightline can be estimated (in units of r^ii), the fitting formulas 
presented in Sect. 5 may be directly applied to assess the mean 
and variance in the expected MACHO optical depth derived from 
a spherical model for the foreground halo. This situation occurs 
when estimating the contribution from non-baryonic MACHOs in 
the dominating halo to the total microlensing optical depth towards 
a gravitational arcs in a galaxy cluster, or an individual image of a 
strongly lensed quasar. As described in Section 5, the error on the 
optical depth due to triaxiality can in this situation easily amount 
to factor of ~ 2 (assuming Eq.[l7]for the optical depth). Note that 
this will be compounded by the surface mass density error coming 
from the uncertainty in the mass and concentration parameter of the 
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spherical model for the foreground halo. In situations where many 
light sources are monitored, and these are projected across a large 
area of the foreground halo (as in the cluster-quasar and cluster- 
cluster microlensing monitoring programmes of Totani 2003 and 
iTadros, Warren & HewetliriOOl , respectively), care must be exer- 
cised when using the formulae presented in Sect. 5. When aver- 
aging over different, but not independent sightlines, the effects of 
triaxiality may be substantially diminished, and the optical depth 
error derived from these formulae should be considered a conser- 
vative upper limit. 

In many cases, however, the impact parameter of the dominat- 
ing halo is not known, and there may be more than one halo giving 
significant contributions to the optical depth of a given sightline. 
This happens when light sources, especially at high redshift, are 
randomly selected without reference to any foreground object. In 
this case, the average microlensing optical depth is often computed 
using: 

_ ^ 3HonMACHO {l + zfDi,D,idz 

where Zs is the redshift of the light source studied and r^MACHO 
is the cosmological density of MACHOs relative to critical at 
zero r edshift. This estimate, often referred to as the lPress & GunnI 
( Il97 approximation, also assumes a constant comoving num- 
ber density of MACHOs, i.e. that the MACHO population does 
not evolve as a function of redshift. Since the matter of the Uni- 
verse is clustered, one does however expect a certain scatter around 
this average, since some sigh tlines will contain more ma tter (and 
hence MACHOs) than others. Zackrisson & Riehml j2007h find, us- 
ing a model that takes into account the clustering of MACHOs into 
spherical CDM haloes and subhaloes, that the distibution of MA- 
CHO optical depths around r is reasonably well described by a log- 
normal function with standard deviation (TinT(2s). As the number 
of intervening haloes increases when more and more distant light 
sources are considered, the sightline-to-sightline scatter, and hence 
o"inr(%), decreases with increasing Zg. In this model, the optical 
depth scatter is dominated by the different number of haloes along 
each sightline, combounded by the different masses, concentration 
parameters and impact parameters for each of these objects. 

Since non-sphericity introduces additional optical depth scat- 
ter on top of that produced by the spherical halos, one may expect 
halo triaxiality to give an significant contribution to (J\i^r{zs), but 
wh en implementing the Scdm (?') scatter formulae derived here in 
the lZackrisson & Riehml j2007h microlensing code, we find the im- 
pact of non-sphericitiy to be negligible. The reason for this is that 
the amplitude of the scatter stemming from triaxiality is relatively 
small compared to that coming from other aspects of CDM clus- 
tering. As can be seen in Figure (2] the mean surface mass density 
varies by more than a factor of ~ 1000 between the innermost 
regions of the halo and rvir, whereas triaxiality generates varia- 
tions of less than a factor of ~ 10 at each impact parameter. This 
means that, once random halo impact parameters are considered, 
the resulting MACHO optical depth distribution will be dominated 
by the scatter introduced by the form of the surface mass density 
profile, while halo triaxiality will be responsible for only a very 
slight modification of this distribution. Considering the fact that 
most sightlines towards high-redshift sources pass inside rvir of 
more than one halo (see Zackrisson & Riehm 2007, for an estimate 
of how many), and that these are likely to have different masses and 
concentration parameters, the impact of triaxiality on the optical 
depth distribution quickly becomes negligibly small. Hence, it can 
safely be ignored in this situation. This greatly reduces the com- 



putational complexity of MACHO microlensing models for high- 
redshift sources. 



7 DISCUSSION 

The results presented here do suffer from a number of shortcomings 
which should be pointed out. The simulations used are dissipation- 
less. In reality, dark matter haloes contain baryons, and the dissi- 
pation and feedback associated with these will inevitably affect the 
overall potential of the system, and thereby the spatial distribution 
of the CDM. According to current models, baryonic cooling will 
increase the central density of the CDM ( e.g. lOnedin et alj|2004) 
and also make the halo more spherical tKazantzidis et al.ll2004) 
The significance of these effects are, however, still difficult to pre- 
dict reliably, as the gas dynamical simulations involved still suffer 
from so-called "overmergin g" problems (e.g. Balogh et al. 200l|; 
ISpringel & Hemauistll2002h . 

When calculating the surface mass density profiles, we have 
moreover considered only the matter present within rvir of each 
halo, whereas simulations have show n that galaxy siz ed CDM 
haloes extend at least out to 2-3rv ir i |PradaetalJ[2003) . We re- 
stricted our analysis to rvir because it becomes increasingly de- 
manding to separate halo particles from the background the further 
one wants to extend the analysis. Even in the presented case, we 
need to separate particles out to ~ l.Srvir because the smoothed 
particles extend their influence inside the virial radius region even 
though they are positioned outside it. We tested our method out to 
3rvir, but the number counts of the halo particles at those distances 
are too low to produce reliable results. Our halo sample does not 
have the resolution needed for extending the analysis further than 
rvir safely. 

The most significant limitation of this paper is the small num- 
ber of haloes in our analysis. This is of course due to the limited 
resolution of the cosmological simulations we had access to. We 
would like to repeat our analysis with a more complete statistical 
sample of haloes, which would hopefully confirm our analytical 
description with smaller error bars. 

We also note that our analytical description of the surface mass 
density is more reliable in the case in which subhaloes are excluded. 
This is because subhaloes can introduce significant mass peaks to 
some radial bins. These peaks can lead to unwanted effects in the 
log-normal fitting procedure which is designed to handle relatively 
smooth and continuous mass distributions within a bin. Large sub- 
halos can also disturb the NFW fits, at last in the low density re- 
gions. Thus, the use of the models which include subhaloes is dis- 
couraged. 



8 SUMMARY 

We have compared a sample of CDM N-body haloes to the smooth, 
spherically symmetric NFW density profile model in three dimen- 
sions. The differences in surface mass density of the haloes and the 
model are studied, and an analytical description of the differences is 
constructed. This description can be used to estimate or reproduce 
the differences between CDM N-body haloes and, in principle, any 
analytical halo model. It can be used in applications in which the 
line-of-sight surface mass densities of CDM haloes play an impor- 
tant role, such as microlensing. 

Our halo sample consists of 27 independent CDM haloes at 
~ 10 redshift snapshots between 0.0 < z < 1.5. The haloes are 



The surface mass densities of CDM haloes 9 



extracted from six cosmological simulations with comoving box 
sizes of 10 h'^Mpc, 40 h~^Mpc and 64 h~^Mpc. The haloes are 
treated both with and without their subhaloes, and the halo sample 
is divided in two mass classes, separated by a mass gap at AI ~ 
10^'^ MqH^^. The analytical description is given for all four cases. 

We find that the surface mass density of the haloes can deviate 
from the spherical model considerably. At minimum, with zero im- 
pact parameter and redshift, the l-cr limits around the NFW surface 
mass density are close to cr = ±20% or a = ±30%, depending 
which haloes are under investigation. At maximum, with impact 
parameter close to r^ii and z = 1.5, the values can be as high as 
a+ = +70% and a- = —80%. The geometric mean of the sur- 
face mass density is offset from the NFW predicted value by —3 % 
to —44 %, depending on the case. 

We also find that the departure from the NFW profile is log- 
normally distributed around the model value. In most cases, the 
median of the surface mass density of the haloes is slightly lower 
than predicted by the NFW profile. The variation of the surface 
mass density around the NFW value grows with increasing impact 
parameter and redshift. 

As an application, we introduce our analytical description to 
the optical depth calculations of MACHOs. In this case, we find 
that the variance in surface mass density due to halo shapes can be 
overwhelmed by the variance caused by random impact parameters 
between halos on the same sightline. 
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